TK!2
=v
=u
=l
#1
:n
astart
:c
ordered list of row-lists of LU matrix  a
:v
'a1
'a2
'a3
'a4
'a5
'a6
'a7
'a8
'a9
'a10
E
#2
:n
a0
:c
list of row-lists of input matrix  a0
:v
'a01
'a02
'a03
'a04
'a05
'a06
'a07
'a08
'a09
'a010
E
#3
:n
a01
:c
row-list
:v
74
-25
66
-11
-96
22
-40
-2
-22
-3
E
#4
:n
a02
:c
   "
:v
58
54
-50
36
-21
-34
23
-64
44
-58
E
#5
:n
a03
:c
   "
:v
-37
-100
-41
-59
-99
52
75
42
99
87
E
#6
:n
a04
:c
   "
:v
81
-68
79
87
-39
55
-62
-45
-2
18
E
#7
:n
a05
:c
   "
:v
-86
17
-58
-71
26
-38
35
-53
90
33
E
#8
:n
a06
:c
   "
:v
-25
81
-93
-78
-11
48
57
80
14
-90
E
#9
:n
a07
:c
   "
:v
-56
95
5
43
68
40
43
-34
-89
20
E
#10
:n
a08
:c
   "
:v
86
46
22
8
25
42
-87
27
-35
99
E
#11
:n
a09
:c
   "
:v
-44
-79
-79
-94
83
68
58
-76
-67
-84
E
#12
:n
a010
:c
   "
:v
-84
-24
56
25
63
22
-55
-24
47
-26
E
#13
:n
A
:c
list of row-lists of the inverse matrix
:v
'A1
'A2
'A3
'A4
'A5
'A6
'A7
'A8
'A9
'A10
E
#14
:n
A1
:c
row-list
:v
E
#15
:n
A2
:c
   "
:v
E
#16
:n
A3
:c
   "
:v
E
#17
:n
A4
:c
   "
:v
E
#18
:n
A5
:c
   "
:v
E
#19
:n
A6
:c
   "
:v
E
#20
:n
A7
:c
   "
:v
E
#21
:n
A8
:c
   "
:v
E
#22
:n
A9
:c
   "
:v
E
#23
:n
A10
:c
   "
:v
E
#24
:n
a
:c
reshuffled list of row-lists of LU matrix  a
:v
'a10
'a7
'a9
'a2
'a3
'a4
'a6
'a1
'a8
'a5
E
#25
:n
a1
:v
E
#26
:n
a2
:v
E
#27
:n
a3
:v
E
#28
:n
a4
:v
E
#29
:n
a5
:v
E
#30
:n
a6
:v
E
#31
:n
a7
:v
E
#32
:n
a8
:v
E
#33
:n
a9
:v
E
#34
:n
a10
:v
E
#35
:n
u
:v
'u1
'u2
'u3
'u4
'u5
'u6
'u7
'u8
'u9
'u10
E
#36
:n
u1
:v
1
2
3
E
#37
:n
u2
:v
0
1
2
E
#38
:n
u3
:v
-1
0
1
E
#39
:n
u4
:v
-2
-1
0
E
#40
:n
u5
:v
E
#41
:n
u6
:v
E
#42
:n
u7
:v
E
#43
:n
u8
:v
E
#44
:n
u9
:v
E
#45
:n
u10
:v
E
#46
:n
v
:v
'v1
'v2
'v3
'v4
'v5
'v6
'v7
'v8
'v9
'v10
E
#47
:n
v1
:v
0
-1
-2
-3
E
#48
:n
v2
:v
1
0
-1
-2
E
#49
:n
v3
:v
2
1
0
-1
E
#50
:n
v4
:v
E
#51
:n
v5
:v
E
#52
:n
v6
:v
E
#53
:n
v7
:v
E
#54
:n
v8
:v
E
#55
:n
v9
:v
E
#56
:n
v10
:v
E
#57
:n
w
:v
'w1
'w2
'w3
'w4
'w5
'w6
'w7
'w8
'w9
'w10
E
#58
:n
w1
:v
E
#59
:n
w2
:v
E
#60
:n
w3
:v
E
#61
:n
w4
:v
E
#62
:n
w5
:v
E
#63
:n
w6
:v
E
#64
:n
w7
:v
E
#65
:n
w8
:v
E
#66
:n
w9
:v
E
#67
:n
w10
:v
E
#68
:n
aux
:v
'aux01
'aux02
'aux03
'aux04
'aux05
'aux06
'aux07
'aux08
'aux09
'aux10
E
=f
#1
:n
LUDMIN
:c
Matrix inversion, matrix decomposition
:t
Procedure
:i
n,a0,a,A
#2
:n
LUD
:c
LU decomposition
:t
Procedure
:i
n,a0,a
:o
D
#3
:n
LUFBSm
:c
Forward and backward substitution for LU matrix inversion
:t
Procedure
:i
n,a,A
#4
:n
Mblank
:c
Blank matrix (i.e. subordinate lists)
:t
Procedure
:i
M,n
#5
:n
Mdelete
:c
Delete matrix (i.e. subordinate lists)
:t
Procedure
:i
M
#6
:n
Mcopy
:c
Copying an  n x n  matrix
:t
Procedure
:i
a0,a
#7
:n
Mproduct
:c
Multiplication of two matrices
:t
Procedure
:i
A,B,C
#8
:n
Transpose
:c
Transposing table A into table B
:t
Procedure
:i
A,B
#1
:n
LUDMIN
:z
" This procedure calls 'LUD' and 'LUFBSm' procedures necessary for inverting
" the matrix a0, and then deletes the auxiliary lists (cancelled here so that
" the LU matrix can be observed).

call LUD(n,a0,a;D)        " LU decomposition
call LUFBSm(n,a,A)        " forward and backward substitution

call Mdelete('a)
call delete('a)
call delete('idx)
/E
7
8
E
#2
:n
LUD
:z
" Notation:  n      - number of equations
"            a0     - lists with given matrix coefficients
"            a      - lists with decomposed matrix coefficients
"            scale  - list with scaling coefficients
"            idx    - list with row pointers (needed only for solving eq-ns)
"            D      - determinant (needed only for determinant evaluation)
"            astart - correctly ordered names of row-lists of  a

" Supporting functions:  Mcopy

" Description: The procedure evaluates the coefficients of lower and upper
" triangular matrices the product of which equals the original matrix  a0 .
" The resulting matrices are combined in matrix  a  (with the unity diagonal
" coefficients of the lower matrix hidden).  The rows of the resulting
" matrix are reshuffled, and one has to refer to list  a  to find their 
" proper order.

" The function 'LUD' may be used for solving linear equations, matrix
" inversion, and determinant evaluation.

" Note: The matrix coefficients have to be stored in lists row by row. 
" The lists  a0  and  a  contain the names of the appropriate row-lists.

call listcopy('astart,a)  " copy ordered listnames into list  a
call Mcopy(a0,a)          " copy matrix coefficients into work matrix  a
zero:= n*1e-32            " floating point zero threshold
D:= 1                     " seed for determinant

for i=1 to n              " generate scale coefficients for pivoting
   big:= max(max(a[i]),-min(a[i]))
   if big<zero then error:= matrix_with_zero_row
   'scale[i]:= 1/big
next i

for j=1 to n                 " column loop of Crout's algorithm
   for i=2 to j-1            " solve for upper triangle coefficients
      sum:= a[i][j]
      for k=1 to i-1
         sum:= sum - a[i][k]*a[k][j]
      next k
      a[i][j]:= sum
   next i
   big:= 0
   for i=j to n              " solve for lower triangle and
      sum:= a[i][j]          " diagonal coefficients
        for k=1 to j-1
           sum:= sum - a[i][k]*a[k][j]
        next k
        a[i][j]:= sum
      pivot:= 'scale[i]*abs(sum)                  " evaluate pivot and com-
      if pivot>=big then (imax,big):= (i,pivot)   " pare it with the last
   next i
   if j=imax goto A
     (a[imax],a[j]):= (a[j],a[imax])       " swap rows
     'scale[imax]:= 'scale[j]              " replace scale
     D:= -D
   A: 'idx[j]:= imax
   if abs(a[j][j])<zero then error:= singular_matrix
   coeff:= 1/a[j][j]
   for i=j+1 to n                    " divide lower triangle coefficients
      a[i][j]:= a[i][j]*coeff        " by pivot elements
   next i
next j
for i=1 to n                         " evaluate determinant
   D:= D*a[i][i]
next i

call delete('scale)
/E
E
#3
:n
LUFBSm
:z
" Notation:  n    - number of rows and columns in the matrix
"            a    - list of row-lists of LU matrix
"            A    - list of row lists of inverse matrix (target)
"            x    - columns of identity matrix converted into columns of
"                   the inverse matrix  A 

" Description: This procedure takes the LU matrix  a , and creates the
" inverse matrix to the original matrix  a0  by column-by-column forward and
" backward substitution using columns  x  of an identity matrix as the right
" hand side coefficients.

for j=1 to n                         " 0-column for identity matrix
   'x0[j]:= 0
next j
for jj=1 to n                        " composing inverse matrix by columns
   call listcopy('x0,'x)
   'x[jj]:= 1
   ii:= 0
   for i=1 to n                      " forward substitution
      k:= 'idx[i]
      sum:= 'x[k]
      'x[k]:= 'x[i]
      if ii=0 goto A
      for j=ii to i-1
         sum:= sum - a[i][j]*'x[j]
      next j
      A: if and(ii=0,sum<>0) then ii:=i
      'x[i]:= sum
   next i
   for i=n to 1 step -1              " backward substitution
      sum:= 'x[i]
      for j=i+1 to n
         sum:= sum - a[i][j]*'x[j]
      next j
      'x[i]:= sum/a[i][i]
      A[i][jj]:= 'x[i]
   next i
next jj
call delete('x0)                     " deleting auxiliary lists
call delete('x)
/E
E
#4
:n
Mblank
:z
" The procedure blanks (i.e. erases the contents) of the first  n  lists
" named in the master list  M .

for i=1 to n
   call blank(M[i])    " error indicates no listname in ith element of M
next i
/E
E
#5
:n
Mdelete
:z
" The procedure deletes subordinate lists of the master list  M  from the
" List Sheet. The master list remains intact.

for i=1 to length(M)
   call delete(M[i])    " error indicates no listname in ith element of M
next i
/E
E
#6
:n
Mcopy
:z
for i=1 to length(a0[1])
   call listcopy(a0[i],a[i])
next i
/E
E
#7
:n
Mproduct
:z
" The procedure generates a product of two matrices defined by the master 
" lists  A  and  B  and their subordinate lists.  The resulting matrix is
" stored in lists subordinate to the master list C, which must be set up in
" advance. Also, the master list  'aux  for temporary storage of a transpose
" of  B  must be available.  The master lists  C  and  'aux  must be at least
" as long as the list B[1]. 

" The matrices  A  and  B  must be compatible, i.e. if  A  represents  m x n
" matrix, B must be  n x m  and the resulting matrix is  n x n .  The pro-
" cedure does not check explicitly the compatibility.  The length of the
" list B[1] (i.e. the 1st row of the  B  matrix) is taken as  n .  Any blanks
" in the lists constituting the matrices  A  and  B  are interpreted as zero
" coefficients by the TK built-in 'dot' function which is used for performing
" the actual calculations.

" Library functions used: Transpose, Mblank, Mdelete.

n:= length(B[1])
call Mblank(C,n)    " error indicates missing elements in list of lists C
call Transpose(B,'aux)
for i=1 to n
   for j=1 to n
      C[i][j]:= dot(A[i],'aux[j])
   next j
next i
call Mdelete('aux)
/E
E
#8
:n
Transpose
:z
" The procedure generates a transpose of a matrix defined by lists subordi- 
" nate to the master list  A .  The list  B , e.g. 'z, must be a master list
" containing the names of subordinate lists, e.g. 'z01, 'z02 and so on.

" Other library function used: Mblank

call Mblank(B,length(A[1]))        " blank all B-lists
for i=1 to length(A)
   a:= A[i]
   for j=1 to length(a)
      call listcopy(a,B[j],j,j,i)
   next j
next i
/E
E
=r
#1
:r
"  ***  Matrix inversion,  LU  decomposition  ***
:s
C
#3
:r
call LUDMIN(10,'a0,'a,'A)
#5
:r
" This model inverts a  10 x 10  matrix using LU decomposition and
:s
C
#6
:r
" subsequent forward and backward substitution.  The build-up of the
:s
C
#7
:r
" inverse matrix  A  may be observed in the interactive table 'invmatrix'
:s
C
#8
:r
" (although switching off the "Display Intermediate Values" option in the
:s
C
#9
:r
" Global Sheet speeds up the solution substantially).  The table 'invmatrix'
:s
C
#10
:r
" may be cleared before the solution, if needed, by a direct call to the
:s
C
#11
:r
" procedure 'Mblank' by typing  / E E  and entering  call Mblank('A,10) .
:s
C
#13
:r
" The product of the original matrix  a0  and the inverted matrix  A  may be
:s
C
#14
:r
" generated by a direct call to the procedure Mproduct by typing  / E E
:s
C
#15
:r
" and entering call Mproduct('a0,'A,'w).  The result may be seen in the
:s
C
#16
:r
" interactive table 'check'; it should be very close to an identity matrix
:s
C
#17
:r
" with  1's  along the diagonal and  0's  everywhere else.
:s
C
=p
=t
#1
:n
matrix
:s
Screen
:v
Horizontal
:f
1
:r
""
:c
""
:t
"Original Matrix"
:u
#1
:l
a01
:w
6
:h
""
#2
:l
a02
:w
6
:h
""
#3
:l
a03
:w
6
:h
""
#4
:l
a04
:w
6
:h
""
#5
:l
a05
:w
6
:h
""
#6
:l
a06
:w
6
:h
""
#7
:l
a07
:w
6
:h
""
#8
:l
a08
:w
6
:h
""
#9
:l
a09
:w
6
:h
""
#10
:l
a010
:w
6
:h
""
E
#2
:n
auxmatrix
:s
Screen
:v
Horizontal
:f
1
:r
""
:c
""
:t
"LU Matrix"
:u
#1
:l
a1
:w
6
:h
""
#2
:l
a2
:w
6
:h
""
#3
:l
a3
:w
6
:h
""
#4
:l
a4
:w
6
:h
""
#5
:l
a5
:w
6
:h
""
#6
:l
a6
:w
6
:h
""
#7
:l
a7
:w
6
:h
""
#8
:l
a8
:w
6
:h
""
#9
:l
a9
:w
6
:h
""
#10
:l
a10
:w
6
:h
""
#11
:w
6
:h
""
E
#3
:n
invmatrix
:s
Screen
:v
Horizontal
:f
1
:r
""
:c
""
:t
"Inverse Matrix"
:u
#1
:l
A1
:w
6
:h
""
#2
:l
A2
:w
6
:h
""
#3
:l
A3
:w
6
:h
""
#4
:l
A4
:w
6
:h
""
#5
:l
A5
:w
6
:h
""
#6
:l
A6
:w
6
:h
""
#7
:l
A7
:w
6
:h
""
#8
:l
A8
:w
6
:h
""
#9
:l
A9
:w
6
:h
""
#10
:l
A10
:w
6
:h
""
E
#4
:n
check
:s
Screen
:v
Horizontal
:f
1
:r
""
:c
""
:t
"Product of Original and Inverse Matrix"
:u
#1
:l
w1
:w
6
:h
""
#2
:l
w2
:w
6
:h
""
#3
:l
w3
:w
6
:h
""
#4
:l
w4
:w
6
:h
""
#5
:l
w5
:w
6
:h
""
#6
:l
w6
:w
6
:h
""
#7
:l
w7
:w
6
:h
""
#8
:l
w8
:w
6
:h
""
#9
:l
w9
:w
6
:h
""
#10
:l
w10
:w
6
:h
""
E
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
